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Abstract. 

A new reactive force field to describe proton diffusion within the solid-oxide fuel 
cell material BaZrOa has been derived. Using a quantum mechanical potential energy 
surface, the parameters of an interatomic potential model to describe hydroxyl groups 
within both pure and yttrium-doped BaZrOu have been determined. Reactivity is then 
incorporated through the use of the empirical valence bond model. Molecular dynamics 
simulations (EVB-MD) have been performed to explore the diffusion of hydrogen using 
a stochastic thermostat and barostat whose equations are extended to the isostress- 
isothcrmal ensemble. In the low concentration limit, the presence of yttrium is found 
not to significantly influence the diffusivity of hydrogen, despite the proton having a 
longer residence time at oxygen adjacent to the dopant. This lack of influence is due 
to the fact that trapping occurs infrequently, even when the proton diffuses through 
octahedra adjacent to the dopant. The activation energy for diffusion is found to be 
0.42 eV, in good agreement with experimental values, though the prefactor is slightly 
underestimated. 
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1. Introduction 

Solid Oxide Fuel Cells (SOFCs) are an important technology for energy production, 
including as part of a hydrogen economy. They are of interest for both oxygen and 
hydrogen ion conduction, where the latter is of particular relevance to trying to achieve 
a reduction in operating temperatures [U [2]. Perovskite structured materials have 
attracted considerable interest for proton conduction, and one of the best candidates is 
BaZrOs- When doped with yttrium on the zirconium B-site, this system achieves one of 
highest levels of proton conductivity known, combined with other favourable properties, 
such as good chemical stability [3J. 

Computational methods offer the opportunity to obtain insights into the atomistic 
details of proton conduction within materials. There have been numerous theoretical 
studies to date of alkaline earth metal zirconate-based perovskites, especially those of 
calcium and barium [U |5]. These two systems provide an interesting contrast since 
while barium zirconate is a cubic perovskite, the calcium analogue exhibits a distorted 
orthorhombic structure. Both cases have been examined quantum mechanically using 
density functional theory [HI El IE]- Here barium zirconate is found not to be cubic 
as it possesses a phonon instability along the R-M edge of the Brillouin zone [9]. This 
represents a failure of present day exchange-correlation functionals within both the Local 
Density and Generalized Gradient Approximations. Although some density functional 
studies have reported no distortion from the cubic structure [10], this has been shown 
to be a result of the lack of semi-core states in the pseudopotential used to describe 
barium [9]. For both CaZrOs and BaZrOs, the doping and proton conduction pathways 
have also been studied quantum mechanically. In the former case, this represents a 
challenging undertaking due to the large number of possible symmetry inequivalent 
proton hops [6]. For BaZrOa, the situation is much simpler with proton migration being 
determined by three motions; two hydroxyl rotations, one about the B-O-B axis and 
one orthogonal to this, where B represents an octahedral cation, and a jump between 
two neighbouring oxygens of the same B-site cation along the edge of the octahedron. 

While the study of proton migration in bulk perovskite materials, including those 
with a single dopant, is readily amenable to quantum mechanical examination, the full 
mapping of the potential energy landscape becomes challenging as defects are introduced 
that lower the symmetry, such as grain boundaries, dislocations, and high concentrations 
of dopants. Hence, there is a need to explore other techniques to examine proton 
diffusion on larger length scales and for more complex configurations. The issue of length 
and timescale can be addressed in part by the use of kinetic Monte Carlo using quantum 
mechanically determined rate constants, as has been performed for CaZrOs [H], or via 
vertex coding, as has been applied to BaZrOa [12j . However, the determination of the 
rate constants can be time consuming and assumptions have to be made regarding the 
transferability of values between distinct environments. 

In contrast to quantum mechanical methods, the use of interatomic potentials 
is readily applicable to large and complex atomistic models of materials. There has 
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already been extensive application of simulation to solid oxide fuel cells using force 
fields based on a shell model description of the polarisable oxide anion. Here both 
bulk and defect properties have been investigated, and the thermodynamics of oxygen 
vacancies, water incorporation and proton trapping by dopants are all accessible. While 
the activation energies for oxygen diffusion within perovskites can also be calculated 
from such force fields, the case of proton diffusion is more problematic. Many of the 
force field simulations of hydrogen incorporation into oxides to date have employed 
the model for the hydroxyl group developed by Saul et al [13] in which this entity 
is described as a covalent molecule, i.e. the oxygen and hydrogen are connected by 
a Coulomb-subtracted Morse potential. Because of the form of this potential, the 
0-H bond cannot be smoothly dissociated and therefore the description is limited 
to equilibrium configurations involving hydroxyl groups. Although there have been 
attempts to describe protons using a formal charge model with a Buckingham repulsion 
[14j . as used widely for metal cations, this approach has only met with partial success 
and does not lead to a broadly applicable model for hydrogen in oxides. 

Recently, van Duin et al [15] derived a parameterisation of the reactive force field 
method, ReaxFF [TB], designed for the simulation of proton diffusion in BaZrOa. Here 
the interactions are described as a sum of many different contributions representing 
terms from covalent bond order, through to van der Waal's and electrostatic energies, 
based on structure dependent charges. The many parameters were determined by 
fitting against a database of quantum mechanical potential energy surfaces spanning the 
component metals, alloys, binary oxides and barium zirconate itself. The resulting model 
was demonstrated to accurately reproduce the experimental activation energy for proton 
diffusion during reactive molecular dynamics [TB] . However, this parameterisation suffers 
from a notable limitation in that BaZrOs is predicted to be highly non-cubic as there are 
six imaginary modes in the phonon spectrum at the F-point [17J. If minimised without 
symmetry constraints, and from a perturbed initial geometry, a stable configuration is 
found in which the unit cell is slightly distorted from cubic to triclinic. More seriously, 
the coordination geometry about zirconium is altered to local symmetry with O- 
Zr-0 bond angles ranging between 74 and 108° instead of all being equal to 90°. As 
a result, many of the distortions observed during proton doping may be the result of 
symmetry breaking in the unstable bulk configuration, rather than as a result of the 
defect itself. The consequences for the activation energy of diffusion and pathways are 
more difficult to assess. 

Although it would be possible to refit the parameterisation of the ReaxFF 
methodology to ensure a cubic ground state for BaZrOs, this is a non-trivial undertaking 
since many of the parameters are strongly correlated. Hence, in the present work our 
objective is to explore an alternative approach to deriving a reactive force field for 
proton diffusion in BaZrOs that is simpler to parameterise. Here we first construct 
an unreactive force field model that reproduces a combination of quantum mechanical 
and experimental data for hydrogen in Y-doped barium zirconate. Reactivity is then 
included through the use of Empirical Valence Bond theory [18] , to couple together the 
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possible states of the proton within the material, as has been successfully employed to 
simulate proton transfer in aqueous systems [19] . After describing the derivation of the 
present model, we present results to illustrate its application to proton diffusion in both 
pure and Y-doped BaZrOa. 

2. Methodology 

2.1. Quantum mechanical calculations 

All quantum mechanical calculations have been performed within the framework of 
Kohn-Sham density functional theory (DFT) using the SIESTA methodology [20] and 
code for linear-scaling Hamiltonian construction. Here the core electrons and nuclei 
are represented through the use of non-local norm- conserving pseudopotentials of the 
modified TrouUier-Martins form [21]. Small core pseudopotentials have been utilised for 
Ba, Zr and Y that were generated for valence electron configurations of 5s^5p^, As^Ap^Ad"^ 
and As'^Ap^Ad^, respectively (i.e. all metal species are dications). Kohn-Sham states for 
the valence electrons are expanded as a linear combination of Pseudo Atomic Orbitals. 
The first basis function for each angular momentum channel is the numerical solution for 
the isolated atom on a logarithmic grid, subject to radial confinement. Soft-confinement 
is employed, according to the method of Junquera et al |22j, with a potential of 100 
Ry and an inner radius equal to 0.95 of the outer boundary. The radius of confinement 
for each orbital is set individually based on the spatial extent of the function. Values 
of this cut-off for bound states vary between 6.0 Bohr for the semi-core states to 10.0 
Bohr for the relatively diffuse 5p levels of Y and Zr. To increase the variational freedom, 
additional radial functions were generated using the split-norm construct, with a norm 
of 0.5 for hydrogen and 0.15 for the remaining atoms, as well as including polarisation 
functions for all atoms (2p, 3d and 4f, for H, O and Ba/Zr/Y, respectively). The cationic 
metals were described using a basis set of double- zeta polarised quality (DZP), while 
oxygen and hydrogen were assigned a triple- zeta, doubly polarised (TZ2P) basis. 

In the SIESTA methodology the electron density is expanded in a uniform grid of 
points in Cartesian space. A kinetic energy cut-off of 600 Ry was used to ensure a high 
degree of convergence with respect to this auxiliary basis set. The electronic structure 
was integrated across the Brillouin zone using a Monkhorst-Pack mesh of dimensions 
controlled by the equivalent real-space cut-off of 12 A according to the scheme of Moreno 
and Soler ^]. Geometry optimisations were performed until atomic forces were less 
than 0.01 eV/A, while the cell stresses were considered converged when below 200 bar, 
though in practice the residual stresses are typically an order of magnitude less than 
this by the time the internal force tolerance is satisfied. In calculations where a net 
charge is present in the unit cell, a neutralising background is applied. 

Previous DFT calculations for BaZrOa have explored a range of exchange- 
correlation functionals, including within the localised basis set approach used in the 
present work [9]. All functionals are found to give similar qualitative results in that 
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this material is predicted to be unstable in the cubic form. The main variation is in 
the magnitude of the lattice parameter obtained. For the current work we employ 
one of a number of recent GGA functionals that appear to offer superior results for 
condensed phases. In particular, we choose to use the AMOS functional [21], which has 
been demonstrated to give good results for lattice parameters, the relative energies of 
polymorphs [25j and, importantly when considering proton incorporation, even water 
clusters [20] • 

Given the distorted nature of the equilibrium structure of BaZrOs with all exchange- 
correlation functionals explored to date, this creates an issue with regard to the 
symmetry of the material. For the experimental cubic system there are only a limited 
number of symmetry unique pathways for proton migration. If the ground state DFT 
structure is used then this symmetry is lost and the number of transition states to 
be located increases significantly. Hence, for reasons of both simplicity and to remain 
faithful to the experimental reality, we have chosen to work with the cubic cell. Because 
the phonon instability arises at the zone boundary and not at the F-point, we can 
suppress the effects of the distortion by using supercells based on odd multiples of the 
primitive unit cell. Thus we perform calculations for cells of size 1x1x1, 3x3x3 and 
5x5x5 in the present work. 

2.2. Force field derivation 

All force field derivation has been performed with the program GULP [271 [28] . using 
least-squares refinement against quantum mechanical structures and energy differences. 
Experimental mechanical properties were also included in the fits to bulk BaO and 
BaZrOa to better constrain the curvature of the potential energy surface. The relaxed 
fitting algorithm [22] has been used in which the structure is optimised at every point 
of the fit so that the change in structure can be used as the observable, rather than the 
forces and stresses. 

Although the shell model has been widely used to simulate oxide phases in force 
field studies, here we chose a rigid ion description for Y-doped BaZrOs. The main 
reason for this choice is the desire to be able to perform efficient molecular dynamics 
simulations of proton diffusion and so the incorporation of a shell model would increase 
the cost by approximately an order of magnitude, regardless of whether using an 
adiabatic or fictitious mass Lagrangian algorithm. Furthermore, because BaZrOs is 
a cubic perovskite, the dipolar oxygen polarisability must be kept low otherwise the 
material would undergo a transformation with a rhombohedral distortion. 

Formal charges are adopted for Ba, Zr, Y and O as an oxide ion, while the short- 
range interaction between these species is given by the Buckingham potential: 



For all interactions we have decided to fix the value of the dispersion coefficient, Cjj, 
to be zero. In the case of the cations this is readily justified, with the possible exception 
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of Ba, due to the low polarisability of the species, while for the oxygen-oxygen interaction 
the relatively uniform cohesive attraction in the high symmetry cubic structure could be 
subsumed into other parameters with no loss of quality. Since the use of formal charges 
probably overestimates the binding contribution from the electrostatic interactions, no 
further attractive contribution is required from the dispersion forces. Furthermore, given 
that a fit is being performed to the results of a GGA calculation, there is no long-range 
van der Waal's contribution to the underlying energy surface. 

Fitting of the metal-oxygen interaction parameters was performed sequentially 
to several structures. Firstly, the Ba-0 interaction was fitted to barium oxide alone 
and then held fixed during the fitting of Zr-0 to barium zirconate. Secondly, having 
determined the metal-oxide potentials, the metal-hydroxide potentials were fitted to 
data for an excess proton in pure BaZrOs in order to constrain this interaction in 
the absence of a dopant. Finally, the remaining metal-oxygen potentials involving 
yttrium were determined based on the quantum mechanical results for a 3x3x3 supercell 
containing a single Y/H defect pair in a range of different configurations, as well as for the 
extreme composition of HBaYOa. Here the bond lengths, hydrogen bonding distances 
and relative energies of the configurations for the supercell were selected as observables. 

To describe the hydroxyl group, we adopt a similar approach to previous works 
by using a molecular mechanics description. Here the oxygen and hydrogen possess 
partial charges that are constrained to sum to the formal overall charge of -1, while 
being Coulombically screened from each other. A harmonic potential is used for the 
short-range intramolecular interaction, while a purely repulsive Lcnnard- Jones potential 
is used to capture the intermolecular repulsion between hydrogen and oxygen. The 
charge distribution within the hydroxyl group and intermolecular terms were fitted 
against the quantum mechanical structures and activation energies for hydroxyl group 
rotation within bulk BaZrOs in the absence of the dopant to ensure these two important 
contributions to proton migration are accurately represented, at least within the limits 
of the underlying DFT functional. In order to better reproduce the relative energetics 
of different Y/H defect configurations it was necessary to distinguish between oxide 
ions in the Zr-O-Zr and Y-O-Zr environments, which will be labelled as 02 and 03 
respectively, when interacting with a neighbouring hydroxyl group. A larger Lennard- 
Jones repulsion was found to be necessary between an adjacent hydroxyl group and the 
oxygen coordinated to yttrium (03). 

In order to allow for proton diffusion, reactivity is included using the Empirical 
Valence Bond (EVB) approach. Here two states, 1 and 2, are coupled by a Hamiltonian: 

_ Hii Hi2 
H21 H22 

The on-diagonal matrix elements, Hn and H22: arc simply given by the conventional 
force field energies of the states being coupled, while the ofi^-diagonal elements control 
the mixing between the states. In the present scenario, the two states represent the 
proton being covalently bonded to two different oxygens at the same time in which the 
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force field attributes of oxide and hydroxy! oxygen are interchanged. This approach has 
been widely used for proton transfer, and other, reactions during simulations, especially 
in liquid water, where models for both an excess proton and the hydroxyl anion have 
been developed [HI [301 [3ll |32l [33] . 

In the case of an excess proton in water, it is well documented that the EVB 
scheme must be extended to the multi-state (MSEVB) variant. Due to the nature of 
the Grotthuss mechanism, the proton cannot be considered strongly localised on two 
water molecules alone and so both the first and second solvation shells often need to 
be included explicitly in the Hamiltonian. In the specific case of perovskite fuel cells 
it is highly improbable that a chain of hydroxyl groups will exist such that multi-state 
coupling is required. Firstly, the concentration of dopants, and therefore protons, is 
usually sufficiently low that two hydrogens will rarely be located on nearest neighbour 
oxygens. Secondly, even in a situation where this did occur, the geometric constraints 
are much greater in the solid-state than in the liquid, and so the ability to attain 
a favourable configuration for formation of a coupled chain of hydroxyls is far more 
limited. Hence, even if the case of a single proton coupling to more than two oxygens is 
contemplated, this is unlikely to occur. Furthermore, proton migration between oxygens 
of the same octahedron is found to occur along the edges and avoid the faces where three 
oxygens would become simultaneously coupled to the same hydrogen. 

Within the EVB method it is necessary to define the off-diagonal matrix elements 
that control the strength of coupling between a hydroxyl group and the oxygen to which 
it is most strongly hydrogen bonded. Several functional forms have been proposed for 
the coupling elements in EVB in terms of the reaction coordinate for proton transfer. 
The important properties that they must possess are a smooth decay to zero as the 
distance between the groups tends to infinity and a derivative of zero with respect 
to the reaction coordinate when the 0-II..0 group is symmetric in order to avoid an 
artificial cusp. Similarly, there are numerous choices for the reaction coordinate. Here 
we take the simple choice of a reaction coordinate, Q, that is defined as the difference 
between the two 0-H bond lengths. We then define the coupling matrix element using 
a Gaussian form: 

Hi,j{r) = X^jexpi-Cij.Q^) (2) 

The two parameters in this expression have been determined as follows. Firstly, 
the exponent for decay, (, was selected such that the coupling constant was negligible 
at the equilibrium geometry for a hydroxyl defect in BaZrOs. This allows the fitting 
of the equilibrium structures and the reactive configurations to be separated. Secondly, 
the prefactor for the coupling. A, is just the difference between the force field activation 
energy at the transition state for proton migration and the value calculated quantum 
mechanically. Hence, this quantity was determined for the symmetric transition state 
of a proton migrating in pure BaZrOa. While potentially the coupling matrix elements 
could depend on the local environment of the hydroxyl group, in the spirit of force 
field transferability we make the assumption that the coupling matrix elements are 
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independent of the nature of the surrounding metal cations, though of course the 
reactivity will differ still due to the underlying force field contribution. 

As a final refinement to the present model, we also include a self-energy contribution 
for the situation when the hydroxyl group is coordinated to Y in order to obtain the 
correct relative energy in comparison to this anion being remote to the dopant. While in 
principle this difference could be captured by the interaction between yttrium and the 
hydroxyl group being assigned a different short-ranged potential, this factor is simpler 
to correct for using the self-energy. 



2.3. Molecular dynamics simulations 

Having determined the reactive force field, we are in a position to simulate proton 
diffusion within BaZrOa. Given the low activation energy for migration and elevated 
temperature of operation of fuel cells, it is appropriate to simulate diffusion using 
molecular dynamics. Because the majority of molecular dynamics codes do not allow for 
the Empirical Valence Bond model (EVB-MD), we have written a new program to handle 
this known as ReaxMD. In this code we have implemented a flexible pressure/stress 
control, which allows one to sample configurations from the NST ensemble. This is 
probably not necessary for the present application, since the simulation box is large and 
the system is stiff, but could be crucial in the simulation of softer materials or of smaller 
boxes. Below we describe in detail the employed algorithm. 

The equations of motion that we employ are a modification suitable for flexible cells 
of those derived in |34J. In that work, several approaches were compared and it was 
found that an optimal choice is to combine a standard NPH scheme with a stochastic 
velocity rescaling thermostat [35] so as to sample the NPT ensemble. The stochastic 
velocity rescaling combines the ergodicity of stochastic schemes with the efficiency of 
global thermostats [36] , while also retaining the possibility of defining a conserved energy 
to control the time step discretisation errors. 

Several schemes are available for flexible cell sampling, NSH, see for example 
[37] |3H1 [39]. Here we build our scheme in a fashion similar to that of [39J, but introduce 
corrections so as to control exactly the sampled ensemble. Defining and Pi the 
positions and momenta of the z-th particle, respectively, h as the matrix of the lattice 
vectors, written as columns a = (/in, /?.2i, ^31), b = (/?.i2, ^^22, ^'-32), c = (/iis, /i23, ^33), 
and p as the matrix of the associated momenta, our equations of motion read: 

Pi , P 

rrii W 



P 
h 



next)^ + 2A;BTL] 

on 

W 

Here W is the barostat mass, the mass of i-th particle, V = det [h], the volume of 
the simulation cell, L is a diagonal tensor with {L)aa = 4: — a, fee is the Boltzmann 
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constant, T the temperature. Pint is the internal stress and Ilext is; 

V 



n 



ext 



hho^P 



cxt 



T\-luT 



(4) 



where Pext is the external load and has to be symmetric to avoid any unphysical rotations 
of the cell. The introduction of a reference cell, ho, is needed here to take into account 
the internal elastic energy accumulated in the simulation cell [37| 138] . In order to avoid 
any cell rotations we fix the reference system so that h is upper half triangular and 
build an initial p that is also upper half triangular. This is equivalent to removing the 
rotational degrees of freedom of the system, as discussed in [39] . 

By inspection, one can prove that the equations of motion ([3]), besides the total 
momentum and the position of the centre of mass, also conserve the following quantity; 

Tr[p^p] ^ 



H = K + U + Uei + 



2W 



2A;Br^(4-a)log(/i), 



(5) 



where K is the total kinetic energy of the particles, U is their potential energy, and Uei 
is the elastic energy stored in the system [38] ; 



U, 



el 



VaTl 



(6) 



where Vq is the volume of the reference cell and e is the strain tensor, which can be 
calculated as: 

e = i((ho^)-^h-hho-^-l). (7) 

By a careful calculation one can show that the Jacobian of the equations of motion ^ 
is proportional to (^)ii(/i)22(^)33- This result is different from that of Martyna et al [39] 
since we are not including a thermostat at this stage and we are fully taking into account 
the constraints on the cell shape (upper half triangular) when computing the Jacobian. 
Combining the above information, one can state that ^ produces the ensemble; 

Vnsh{Pi r, h, p)d]9drd^hd^p oc 



niiTi 



X 6 



J2Pi ]^iH- Ho) dpdrd^hd^p 



where Hq is the initial value of H. Here we write d^h as a reminder that we 
are only sampling the 6-dimensional space of all the possible upper half triangular 
cell matrices. The change of measure can be performed by considering the relation 

(/i)?i(/i)22d^hocd9h. 

When combined with a thermostat, the resulting ensemble is: 



VMsrip, r, h, p)dpdrd^hd^p oc VS 



knT 



(9) 



dpdrd^hd^p 
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This ensemble is the usual NST ensemble, with an extra V term that compensates for 
the constraint on the centre of mass of the system. In this way, the expectation value of 
any observable depending only on the inter-particle distance or cell shape is unaffected 
by the centre of mass constraint and is identical to what would be obtained from a 
Monte Carlo simulation. For a more detailed discussion of this term, see [31]. The 
apparently asymmetric contribution in ^ is indeed necessary to explore the correct 
ensemble. We have tested our implementation on a small model system and verified 
that, for an external isotropic pressure, this correction is needed to recover invariance 
with respect to permutation of the lattice vectors. 

In order to numerically forward integrate the equations of motion (jsj) for a time 
interval. At, we use a scheme based on the Trotter-decomposition jlOl HI], which is a 
straightforward extension of the algorithm used for the isotropic, NPT, case [34]: 

• Evolve thermostat equations (T-step) for At/2. 

• Evolve velocities (P-step for At/2. 

• Evolve positions and velocities (i?-step) for At. 

• Evolve velocities (P-step) for At/2. 

• Evolve thermostat equations (T-step) for At/2. 

The T-step rescales the particle velocities and the cell momenta with the same factor, 
which is obtained through a stochastic procedure [3^ 135] . Note that here the total 
kinetic energy now includes also the cell kinetic energy, Tr[p-^p]/2IV, and that the 
number of degrees of freedom is augmented so as to also account for the 6 cell degrees 
of freedom. 

The P-step is the evolution of the atoms' and cell momenta and can be written as; 

At 




{P)ap{t) + [l^(Pint - next) + 2A;bTL 



At 

T 



2mi \ 2 J 

(f-(t))„(f,(t))^ fAtV 



Srrii V 2 / 

where the last two terms on the right hand side represent the time propagation of the 
kinetic contribution to the stress tensor. The instantaneous internal strain. Pint, is 
calculated here as; 



{Pmt)ap(t) - — 



iU'] 



dU{r,h) 



and setting to zero the terms below the diagonal: 

(Pint)a/3 = for a > /3 (12) 
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The Q-step is the evolution of the atoms' positions and cell vectors and can be 
written as: 



p(t) 



Vi{t + At) = ei^^*ri(t) + 



\pit)] 


-1 







p(t) 
e v^' 
p(t) 



p(t) 
e vv 



At 



— e 



p(t) 



At 



nii 



At 



(13) 



P^it + At) 

h{t + At) =ei^^*h(t). 

The matrix e^^* is obtained from ^^At by means of the Crank-Nicolson expansion 
truncated at the 5^^ order |42j . 

When a finite time step is adopted, sampling errors are inevitably introduced. To 
keep them under control, one can monitor the drift of an effective energy that quantifies 
the detailed balance violation [33113]. Also in the present case, the effective energy drift 
is computed by book-keeping the total kinetic energy increment given by the thermostat 
and subtracting it from ([s]). 

All molecular dynamics simulations have been performed with a stochastic 
thermostat and the above new barostat with relaxation times of 0.1 ps and 1 ps, 
respectively. All simulations were 1 ns long, unless otherwise stated, with a 1 fs time 
step that led to perfect energy conservation for non-reactive trajectories, while for the 
reactive ones at high temperature a slight drift was observed ( 3 eV out of a total 
absolute magnitude of 31200 eV, which corresponds to a 0.01% variation, or a drift of 
1 meV per degree of freedom per ns). For every temperature, 12 runs independent runs 
were performed, employing different starting configurations and random velocities. 



3. Results and discussions 

3.1. Quantum mechanical calculations 

Optimisations of both BaO and BaZrOs have been performed using the AMOS exchange- 
correlation functional. As shown in Table [l| the present results show very good 
agreement with the experimental cell parameters, especially for BaZrOs. While 
the many previous studies of this material have typically exhibited the systematic 
overestimation/underestimation of the lattice parameter associated with use of a 
GGA/LDA functional, the AMOS functional appears to be superior for the present 
structures. We note that the inclusion of the f polarisation functions for Ba and Zr are 
important here, in that the lattice parameter increases to 4.210 A from 4.1898 A in the 
absence of these. For BaO, the error in the lattice parameter for AMOS is larger, but 
still only 0.3 %. 

There have been several previous quantum mechanical studies of barium zirconate 
and the results found in the present work are consistent with the more recent of these [9] 
in finding that the structure is phonon stable at the F-point, but exhibits instabilities 
elsewhere in the Brillouin zone when cubic. As a result, we have performed defect 
calculations at constant volume and using odd multiples of the primitive cell in order 
to suppress the influence of these unstable modes. 
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Before considering the doping of BaZrOa, we first contemplate the incorporation of 
a proton in a single unit cell of pure material. Here a charge neutralising background 
is applied to compensate for the positive charge of the proton. There is only one 
symmetry unique structure for the proton in this cell in which the hydroxyl group lies at 
a small angle away from being parallel to a lattice vector and directed partly towards a 
neighbouring oxide ion in the same plane as the common zirconium ion. From this site, 
there are two distinct types of activated motion that can occur in which the hydroxyl 
rotates about two orthogonal axes, one leading to a new intermolecular hydrogen bond 
to an oxygen within the same octahedron and the other leading to interaction with 
a different octahedron. Of these two processes, the first is the most energetic with a 
computed activation energy of 0.246 eV. 

In addition to computing the pathways for hydroxyl group reorientation, it is also 
possible to calculate the barrier to proton transfer between two adjacent oxygens. By 
exploiting the symmetry of the transition state, this configuration can be readily located 
via a constrained minimisation. Here a low barrier of 0.0375 eV is computed, suggesting 
that proton exchange between neighbouring oxygens may be faster than rotational 
reorientation, though we will discuss a number of caveats on this below. 

The activation energies for proton migration in regions away from the dopant have 
also been calculated by Bjorketun et al. [H] using the PW91 functional Here they 
found barriers of 0.16-0.18 and 0.20 eV for rotation and transfer, respectively. Gomez et 
al. [To] have also computed the same quantities, with the same functional, but obtain 
slightly different values of 0.14 and 0.25 eV. Qualitatively both of these sets appear to be 
different to the present work in that the order is reversed, i.e. rotation would be easier 
than proton transfer. One of the key differences between the literature studies and the 
present work is the functional used and thereby the lattice parameter for BaZrOs. Using 
the same computational conditions, we have computed the variation of the rotational 
energy barrier by increasing and decreasing the lattice parameter. For values of 4.0 
and 4.4 A we obtain values of 0.513 and 0.170 eV, respectively, thus illustrating the 
strong dependence on lattice parameter. Consequently, the lower value for this quantity 
found in earlier work is entirely consistent with their larger lattice constant of 4.25 A. 
Conversely, the activation energy for the proton transfer reaction is found to increase 
rapidly with increasing lattice parameter, again in line with the difference between the 
studies. 

A study of proton migration in Y-doped BaZrOs by Merinov and Goddard HI 
[16] also examines the same processes, this time using calculations based on the PBE 
functional [17] . Given that both PBE and PW91 are GGA functionals, one might be 
expecting them to yield similar results for the present material (though we note that 
this is not always the case [IE])- However, a higher activation barrier for proton transfer 
of 0.48 eV is quoted. In this work they also note the strong sensitivity of the barrier to 
the 0-0 distance along the edge of the octahedron, but there is no explicit information 
given regarding the optimised lattice parameter. Thus it is difficult to explain the 
apparent discrepancy between this work and the earlier works in the field based on this. 
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However, a recent study by Gomez et al. [I2] using the PBE functional also obtains 
proton transfer barriers of 0.5-0.6 eV away from the dopant. An even higher activation 
energy for proton transfer of 0.69 eV has been determined for BaZrOs by Miinch et al. 
[1] . However, this value was computed within the Local Density Approximation and so 
systematic overbinding is to be expected. 

Given the self-interaction errors inherent in GGA Kohn-Sham theory, any values for 
proton migration activation energies are likely to be an underestimate, but the higher 
value found in the present work relative to Bjorketun et al. is arguably preferable on 
the basis of the more accurate lattice parameter. The fact that Merinov and Goddard 
III obtain closer agreement with experiment may not be indicative of a better result, 
given the expectation that Kohn-Sham theory with a GGA functional should not be 
perfect. Indeed, without considering the influence of thermal fluctuations on the system 
one might not expect to agree with high temperature experimental data. 

An additional reason for the much lower barrier to migration in the present work, as 
compared to those of Merinov and Goddard HI or Gomez et al. relates to the supercell 
dimensions. In both of the aforementioned studies a 2x2x2 supercell was employed, 
which allows the barium zirconate structure to distort due to the imaginary phonons 
at the zone boundary, whereas in our work we restrict this by using odd numbered 
supercells. Therefore, in the higher literature values there is a contribution to the 
activation energy from reorientation of the octahedra, which is absent in the current 
study, and we presume in the experimental system that is actually cubic. Furthermore, 
the low barrier here is for a single cell, rather than a supercell and so image interactions 
may contribute to lowering the magnitude. While this implies that the present value is 
not an accurate estimate of the barrier for proton migration in the low concentration 
limit, the objective here is only to provide a potential energy surface to train a force 
field against, and so provided the fit is performed for the same configuration this image 
interaction is not an issue. The barrier height for more dilute configurations can then 
be determined subsequently using the force field instead. 

Next we turn to consider the question of proton-dopant binding for Y in BaZrOa. 
Here we have examined five different proton binding configurations, shown schematically 
in Figure [T] In the present work, there are two distinct minima for a proton bound to the 
same oxygen within the plane of each permutation of pairs of lattice vectors (i.e. the OH 
group is not orthogonal to each B1-0-B2 direction, where B1/B2 = Y or Zr, but is tilted 
slightly towards Bl or B2). For this reason, there are a greater number of minima here 
than in the work of Bjorketun et al.. The most stable state is for the proton to be bound 
to an oxygen coordinated to yttrium, with the hydroxyl tilted to interact with another 
oxygen of the Y octahedron. In Table |2] the energies of all configurations relative to this 
state are presented. Earlier DFT calculations suggest that the proton prefers to bind 
to an oxygen between two Zr ions, rather than adjacent to Y, albeit with a negligible 
difference of 0.01 eV. Here we find that same magnitude of energy difference, but with 
the opposite sign. In a calculation where the polarisation functions are omitted on Ba 
and Zr, leading to a large cell parameter, this ordering is again reversed, suggesting that 
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this value is also slightly sensitive to cell volume. Given that the energy differences are 
of order of thermal energies per degree of freedom at ambient conditions, it appears that 
the two sites are effectively isoenergetic. 

Comparing the energy difference between the most stable proton binding site and 
that most remote from yttrium within the supercell provides an estimate for the trapping 
energy of -0.237 eV. Stokes and Islam |19] have computed this quantity based on shell 
model potentials and obtained a value that appears to be similar, but marginally more 
exothermic (value is estimated from Figure 3 in their work). As they note, there are no 
corresponding experimental values for exactly this system, though the value lies within 
the range of -0.2 to -0.4 eV measured for other materials. Although a recent neutron 
spin-echo study [50] has examined the proton dynamics in 10% Y-doped BaZrOs, no 
value for the trapping energy is quoted. 

3.2. Force field 

The quantum mechanical information from the previous section has been used as the 
basis of a force field derivation for Y-doped BaZrOa, as described in the methods section. 
The final parameters are given in Table [3] Formal-charged shell model parameters 
already exist for BaZrOs [IH]. Although derived independently of the work of Stokes 
and Islam, it is reassuring that the resulting Ba-0 potential is very similar, as is the 
predicted lattice parameter for this material, though they were fitted to different values. 
In the present work the repulsive Zr-0 potential is a little steeper, but this may reflect 
the absence of the oxygen-oxygen attractive contribution. When fitting the interaction 
between the metal cations and the hydroxyl oxygen it was generally found to be sufficient 
to only vary the A term, rather than the p value as well. 

The final results for the fitted force field are compared against the quantum 
mechanical data, experiment, and the results of the ReaxFF force field for the bulk 
properties in Table [Tj Unsurprisingly, the AM05 structures of both BaO and BaZrOa are 
perfectly reproduced as they were heavily weighted in the fit. In contrast, the 
ReaxFF parameterisation overestimates the lattice parameters for both systems when 
constrained to be cubic. With the smaller unit cell dimensions and formal charges, 
the present force field overestimates the hardness of BaZrOs, though it accurately 
captures the on-diagonal elastic constant component, Cu. Because of the high 
symmetry cubic structure, even inclusion of a dipolar shell model would not soften the 
structure. The symmetry of the Born effective charge tensor for oxygen indicates that a 
quadrupolar or higher polarisability model would be required to improve the description 
of the curvature. Conversely, the ReaxFF model predicts mechanical properties for 
BaZrOa that are too soft. This is almost certainly due to the lower charges of +1.66 
and +1.79 for Ba and Zr, respectively, generated by the electronegativity equalisation 
procedure. Although the hardness is a potential weakness of the present model, it 
should have minimal impact on the proton conductivity unless the pressure dependence 
is examined. 
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Having determined the force field parameters required to simulate BaZrOs, it is 
now possible to validate our algorithm for the integration of the equations of motion 
(|3| using this specific case. To do this, we compare the elastic constants of BaZrOs as 
computed from the cell fluctuations [38] recorded during a molecular dynamics run at 
constant temperature and pressure against those calculated by lattice dynamics using 
analytical second derivatives, as implemented within the program GULP [28]. The 
elastic constants are determined from the ensemble average of the strain correlations 
using the methodology of Parrinello and Rahman [51] . A set of runs were performed with 
different thermostat and and barostat relaxation times, hereafter indicated as and r^, 
respectively. As a further check we also calculated the same properties using a Langevin 
thermostat with different values of the friction coefficient. We performed 1 ns long NST 
runs at 300 K and note that for ti, > 10 ps the cell oscillations were extremely slow; 
thus the simulation time was insufficient to achieve complete convergence of the elastic 
constants. On the other hand, < 0.1 ps leads to instabilities in the integration of the 
equation, preventing good energy conservation from occurring with a 1 fs timestep, while 
for intermediate values the calculated elastic constants appear to be quite insensitive to 
the thermostat and barostat relaxation times. 

The calculated elastic constants from the molecular dynamics simulations are 
compared against those from lattice dynamics in Table |4} Here the results obtained 
using the stochastic and Langevin thermostats are essentially equivalent to within 
the statistical uncertainty, with the possible exception of C12. In contrast, there is 
a definite discrepancy between the values of Cu obtained using molecular dynamics and 
that obtained from lattice dynamics for the energy minimised unit cell. However, this 
difference is a result of finite temperature effects; as the unit cell of BaZrOs expands, 
so the value of Cu decreases. If the lattice dynamics calculation is performed using 
the average cell dimensions found during the molecular dynamics trajectory at 300 
K, then the agreement between the techniques is significantly improved. We note 
that the thermal expansion is considerably overestimated by molecular dynamics below 
the Debye temperature due to the lack of vibrational quantisation. Hence, if a free 
minimisation is used to determine the thermal expansion, then the reduction in Cu 
would be considerably less. 

Turning now to consider the energetics of proton defects in Y-doped BaZrOs, 
the results from the present force field are compared against the reference quantum 
mechanical data in Table |2} While the long-range limit of the proton trapping energy is 
explicitly fitted via the self-energy of the hydroxyl group when coordinated to yttrium, 
the relative energies of intermediate configurations are also well reproduced. The 
greatest difficulty is in getting the energetic balance right between subtly different 
configurations where the proton is close to yttrium. Although the percentage errors 
are high for configurations 2 and 3, the absolute error in both cases is less than 0.05 
eV, which is of the order of ambient thermal energy per degree of freedom under the 
conditions at which fuel cells typically operate. Importantly, the ordering of the energies 
is correct and the trapping energy matches the quantum mechanical result. Given that 
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the largest thermodynamic difference to be overcome for proton diffusion away from 
the dopant is correct, it is hkely that the small errors in the proton dynamics close to 
yttrium will not prove critical except at dopant concentrations higher than those used 
experimentally. 

For the final step in the force field parameterisation it is necessary to consider the 
proton transfer event between two adjacent oxygens. Here we have fitted the parameters 
of the coupling matrix elements of the EVB for the case of pure BaZrOs. The resulting 
values of the prefactor (A) and exponent {() are 0.7998 eV and 16.0 A~^, respectively. 
With these values the density functional barrier to proton transfer is reproduced. The 
barrier for hydroxyl rotation as given by the force field is not guaranteed to be exactly 
equivalent to the density functional value since the hydroxyl model has to balance many 
constraints. However, the computed value of 0.231 eV is in good agreement with the 
first principles value of 0.246 eV. 

3.3. Proton diffusion 

Having determined a force field representation of the underlying quantum mechanical 
energy surface, it is now possible to perform molecular dynamics simulations of the 
proton within BaZrOs. Before introducing reactivity into the system, the dynamics 
of a proton in undoped BaZrOs have been examined to check the behaviour of the 
present model. As shown in Figure [2| at 1500 K the proton is able to rapidly explore 
each of the 8 isoenergetic minima accessible to it. In Figure [3] the hydrogen location 
probability isosurface is also shown for a low temperature run at 500 K. Here the eight 
lobes associated with the hydrogen being bonded to the central oxygen are broadened 
into four elongated surfaces. This indicates that even at lower temperatures the hopping 
between different tilt angles of the hydroxyl group occurs rapidly, while the rotation 
about the Zr-O-Zr axis is hindered. 

Activating the Empirical Valence Bond model on top of the underlying potential 
energy surface allows us to now explore the diffusivity of hydrogen in BaZrOs. Here we 
consider the proton migration in both pure BaZrOs and in the presence of an yttrium 
dopant. For the present calculations we have used a 6x6x6 supercell contain 216 formula 
units. Note that the force field model has no imaginary modes anywhere in the Brillouin 
zone and therefore the use of an even-sized supercell will have no repercussions. In 
the present work we consider just a single Y atom within this supercell leading to a 
concentration of 0.46%. While experimental studies typically consider much higher 
concentrations, as do other simulations, we defer the consideration of multiple dopants 
and microstructural complexity to future work. For BaZrOs, the approximation of 
examining a single dopant atom can be justified based on the experimental observation 
that the activation for diffusion is relatively insensitive to the yttrium concentration up 
15% [3]. 

In Figure |4| a sample trajectory for proton diffusion at 1500 K is illustrated. 
Although it is difficult to capture the full dynamics in a static projection, the tendency 
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of the hydrogen to diffuse along the edges of octahedra can be seen, as can the high 
degree of mobihty at this temperature. In this particular trajectory the proton was 
initially placed at a site remote from the yttrium dopant, but the hydrogen diffuses 
past this octahedron during the timescale of the simulation. The low density of time- 
averaged hydrogen positions close to yttrium illustrates that the proton is not trapped 
in this region, but is able to diffuse by readily. This is an important feature of the 
present model. In an earlier version of the model, the use of formal charges led to an 
overestimation of the proton-dopant binding and this resulted in irreversible trapping 
on the timescale of the simulations. Through the introduction of a self-energy term and 
careful fitting of the model to the quantum mechanical energies for a range of proton 
binding sites this overbinding is avoided. 

Although the proton is found to be able to diffuse close to yttrium without becoming 
bound, in some trajectories trapping is observed. At 1500 K, where the largest number 
of events are seen and therefore the statistical reliability is greatest, a proton spends on 
average 5 ps coordinated to a given oxygen between two zirconium ions before migrating 
to an adjacent site. In contrast, when the proton does become bound to an oxygen 
in a Y-O-Zr linkage, the proton resides there for between 50 and 90 ps in all cases. 
Therefore, the residence time is an order of magnitude greater adjacent to the dopant 
at this elevated temperature. At lower temperatures the residence times are longer, as 
would be expected, but the longer lifetime of the proton bound adjacent to yttrium 
persists. 

In order to quantify the rate of proton migration, we have used an expression for 
the diffusion coefficient based on the jump frequency between sites; 

D = -a'u (14) 
6 

where a is the hopping distance and u is the average hopping time measured during the 
simulations. Given that the hydrogen is free to rotate around the octahedra vertices, we 
assumed the hopping length to be the average oxygen distance. Here we have performed 
12 independent simulations for each temperature, where each run samples 1 ns of real 
time. The values from the ensemble of simulations are used to determine the average 
diffusion coefficient, as shown in Figure [5} and the error bar reflects the maximum 
variation between different individual runs. As the frequency of site jumps increases 
with temperature, so the variability between runs decreases as the temperature rises, 
leading to a reduced uncertainty. The data for 500 K are not included in the fits 
to determine the activation energies, since only 25% of all simulations exhibited any 
diffusion events. The result for this temperature is shown as an inset in Figure [5] in 
order to demonstrate that the linear fit to the diffusion coefficients does pass within the 
range of the average value for this temperature and the upper bound to the observed 
migration rate. 

In Figure |5| the results for both doped and undoped BaZrOs are given. As can be 
seen, the two curves are parallel to each other, leading to the same activation energy, 
and are essentially equivalent within the uncertainty. This indicates that the presence 



Reactive Force Field for Proton Diffusion in BaZrOs 



18 



of yttrium at low concentration plays little role in influencing proton diffusivity. This 
arises primarily because the proton is infrequently trapped adjacent to yttrium, and 
even when this occurs diffusion is only halted for several tens of picoseconds, out of a 
nanosecond of simulation, before it escapes. The reduced hydrogen mobility at higher 
yttrium concentrations can therefore be postulated to arise from the loss of a pathway 
through the system that avoids frequent encounters with yttrium. 

Finally, we can compare the proton diffusion coefficients obtained in the present 
work with those of other simulations with the ReaxFF force field and also against 
experiment, as shown in Figure [6j Again, we note in comparing results the caveat 
regarding the different yttrium concentration in the present work. This aside, it is found 
that the activation energy for proton migration is very similar in all cases, with the EVB 
method being marginally closer to experiment. Both the ReaxFF and EVB approaches 
were designed to be parametric representations of a predetermined potential energy 
surface computed using density functional theory, though with very different underlying 
functional forms. Hence, it is satisfying that both methods yield comparable results 
as it suggests that the quantum mechanical energy surface is being captured properly 
overall. 

Where there is a discrepancy between theory and experiment is in the prefactor, 
with both molecular dynamics simulations underestimating the experimental value. 
There are several possible reasons for this, including the neglect of quantum effects for 
the hydrogen nucleus and/or failure to properly describe the vibrations in key regions 
of the energy surface. Interestingly, the combined experimental and theoretical study 
of Karlsson et al. [50] includes first principles estimated diffusion coefficients that have 
a prefactor that is too large relative to the measured values. Given the sensitivity 
of the activation energies to lattice parameter, it may be that thermal expansion and 
the consequences for the underlying energy surface is the key to the prefactor. In 
the case of the first principles results the lattice parameter is the static one, whereas 
in the molecular dynamics simulations the unit cell will expand considerably, and in 
all probability exceed the experimental value due to the lack of quantisation of the 
vibrations. Further investigation of the infiuence of the cell fiuctuations, thermal 
expansion and their correlation with the diffusion coefficient is required to determine 
whether this is a significant factor. 

4. Conclusions 

A new reactive force field model for the diffusion of hydrogen in the solid oxide fuel 
cell material, barium zirconate, has been created based on a potential energy surface 
generated using density functional theory data with the AMOS exchange-correlation 
functional. By combining a conventional rigid ion interatomic potential set with the 
empirical valence bond approach, the process of deriving the force field to describe 
equilibrium configurations and transition states for hydroxyl group rotation can be 
separated from the parameterisation of the proton hop between oxygens along the edge 
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of a B-site octahedron. 

By using AMOS as the functional, the present model is better able to reproduce 
the equilibrium structure of BaZrOa than previous quantum mechanically derived force 
fields, and the ground state is fitted to be cubic in accord with experiment [52] within 
the force field. Although the computed activation energies for hydroxyl rotation and 
proton hopping between two oxygens are quite different to the values from some previous 
quantum mechanical studies, the magnitude of the highest barrier is in line with some of 
the literature values. The origin of this overall similarity can, at least partly, be traced 
back to the variation of the activation energy with cell parameter and the fact that the 
barriers for the two key processes, rotation and hopping, change in opposite directions. 
For the present method, changing the cell length from 4.0 to 4.4 A causes the barrier to 
rotation to decrease from 0.51 to 0.17 eV, while the hopping barrier increases from close 
to zero to 0.31 eV. Hence, as the cell varies in size there will always be a barrier greater 
than 0.2 eV, but the character of the rate limiting process will vary. By fitting a model 
that accurately captures the volume and using molecular dynamics to sample the rate 
of diffusion, while allowing the cell to fluctuate in an NPT ensemble, this may achieve 
the correct balance between different proton motions. This is supported by the finding 
that the overall activation energy, as computed from the temperature dependence of the 
diffusion coefficient, is in good agreement with experiment. 

Having validated a new reactive model for proton diffusion, it is now possible to 
consider the broader application. Although only a low concentration of yttrium has 
been examined in simulations so far, the extension to higher dopant levels should offer 
no complications given that the model includes information on the fully doped limit. 
Based on the comparison of proton migration in pure and doped BaZrOa, albeit at a 
low concentration so far, the effect of the presence of yttrium is minimal on the overall 
diffusion rate as the increased residence time of the hydrogen near the dopant is limited 
by the low probability to become trapped at this site. 

Aside from the consideration of dopant configurations as a function of 
concentration, the use of a reactive force field makes it feasible to examine more 
complex environments, including extended defects, such as surfaces, dislocations and 
grain boundaries, though the present model will require more complex extension unless 
the proton is the only reactive species. Here systematic enumeration of transition 
states through the use of direct quantum mechanical calculations can rapidly become 
prohibitive. Examination of the contribution of these features to the macroscopic rate 
of proton diffusion is now possible in the future. 
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Poisson's ratio (GPa) 
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Table 1. Comparison of structure and properties obtained for BaO and BaZrOs as 
measured experimentally and calculated with two force fields, namely the ReaxFF 
model and that of the present work fitted to the quantum mechanical data. In the case 
of the ReaxFF force field values are given for the unstable cubic BaZrOa structure. The 
optimised cell parameters computed using density functional theory with the AM05 
functional are also given. 
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1 
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0.000 


0.000 


2 


3.110 


0.061 
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3.643 


0.006 


0.047 
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5.584 


0.129 


0.124 


5 


8.979 


0.237 


0.238 



Table 2. Energies of the different configurations for proton in Y-doped BaZrOs, taken 
relative to the most stable state. Both the values given by the quantum mechanical 
calculations with the AMOS functional and by the resulting EVB force field fit are 

reported. The minimum yttrium-hydrogen distance is also given for the optimised 
quantum mechanical geometry in a 3x3x3 supercell as an indication of the proton 
location. 
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Table 3. Intermolecular force field parameters for the rigid ion model derived in the 

present work. Here the atom types 01, 02 and 03 distinguish between the oxygens 
of a hydroxyl group, an oxide ion bonded to Zr only, and an oxide ion bond to Y and 
Zr, respectively. Where just O is given this implies that the potential acts between all 
types of oxygen. Formal charges are used for all species, except H and 01, which have 
charges of 0.308698 and -1.308698 a.u., respectively. All intermolecular interactions 
are of the form of the Buckingham potential, except for the intermolecular H-01 
interaction, which is given by a Lennard-Jones 12-6 potential. The intramolecular 
H-01 interaction is represented via a Coulomb-subtracted harmonic potential with a 
force constant and equilibrium distance of 46.016 eV A^^ and 0.985357 A, respectively. 
Interactions are smoothly tapered to zero over a range of 2 A, starting at 8 A. 
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Stochastic 


Langcvin 


Static (OK) 


Static (300K) 


Cu 


(GPa) 


361 ± 15 


356 ± 7 


371.1 


365.7 


C44 


(GPa) 


137 ± 7 


134 ± 3 


135.3 


134.9 


C12 


(GPa) 


131 ± 15 


139 ± 6 


135.3 


134.0 



Table 4. Average clastic calculated from the cell fluctuations during an NST run 
at 300K with r(,=0.1 ps using different Tj and ( for the Stochastic and Langevin 
thermostats. The values obtained using lattice dynamics (static) with GULP are 
also shown for comparison. Static values are calculated using the optimised lattice 
parameters (0 K) and the average cell dimensions from the molecular dynamics (300 
K). 
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Figure 1. Schematic illustration of the five proton configurations examined relative 
to the position of yttrium. The numbered sites 1 to 5 refer to the configurations in 
Table [2] Here the hydro gens in positions 1 and 2 are bonded to 01, positions 3 and 
4 to 02, while position 5 is bound to an oxygen one conventional lattice parameter 
behind 03 into the plane of the figure. 
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Figure 2. (Colour online) Positions of a non-reactive proton (small white spheres) 
sampled during a short simulation at 1500 K. It is evident how the proton can freely 
rotate around the oxygen to which it is bonded. Here Zr (grey) and O (red) atoms 
are in the centre and on the vertexes of the octahedra and Ba (green) atoms are in 
between them. 
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Figure 3. (Colour online) Probability isosurface for a proton (light grey) during an 
EVB-MD simulation at 500 K. The isosurface shows a four-fold symmetry around a 
Zr-O-Zr bond. Ba (green) atoms are also shown as larger light grey spheres without 
bonds. 
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Figure 4. (Colour online) Diffusion trajectory of a proton (small white spheres) at 
1500 K in the BaZrOs lattice. For clarity the proton trajectory has been smoothed 
and every sphere represents the proton position averaged over five successive frames 
taken at 1 ps intervals. The O (red) atoms form square lattice and the Ba (green) 
atoms are in the centre of the squares. From this view point the Zr atoms are screened 
by the oxygen atoms at the squares vertexes. In the centre of the image a light grey 
atom on the square lattice indicates the oxygen that is nearest neighbours of a Y atom 
underneath. 
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Figure 5. Proton diffusion coefficient at different temperatures for the pure and Y 
doped BaZrOa material. For clarity the low temperature point (500K) has been placed 
in the inset. 
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Figure 6. Comparison between the proton diffusion coefficient calculated in 
this work using EVB-MD (red dots) with the experimental values reported for 
BaZro.gYo.iOs-i by Kreuer 3J (blue squares), Karlsson et al [SD] (purple crosses) 
and with previous simulation of BaZro.875Yo.i25Ho.i25 03 with the ReaxFF force field 
[T5] (black triangles). 



